load("replication_data.RData")
library(tidyverse)
library(clubSandwich)
library(MatchIt)
library(ggeffects)
library(car)
library(stargazer)

# GRAPH PARAMETERS
title.size <- 20
nom.label.size <- 10
x.axis.tick.size <- 14
y.axis.tick.size <- 14
x.axis.label.size <- 18
y.axis.label.size <- 18
strip.text.size <- 14
line.width <- 1.5
axis.text <- 1.2
label.text <- 1.2
text.size <- 1.4

#### Figure 1 ####
ggplot(total_cases,
       aes(x = term, y = total_cases_on_docket)) + theme_bw() + 
  geom_line(size = 1) +
  labs(x="Term", y="Number of Cases", title="") +
  scale_x_continuous(breaks = seq(1880, 2020, 10)) +
  theme(plot.title = element_text(size =title.size, face = "bold", hjust=.5), 
        axis.text.x = element_text(size=x.axis.tick.size),  axis.title.x=element_blank(), 
        axis.text.y = element_text(size=y.axis.tick.size),  axis.title.y = element_text(size=y.axis.label.size),
        plot.margin=unit(c(0,3,0,3),"pt"))

#### Figure 2 ####
ggplot(appeals_and_cert ,  aes(x = term, y = total_cert_grants)) + theme_bw() + 
  geom_line(size = 1) +
  labs(x="Term", y="Number of Cert Petitions Granted", title="") +
  scale_x_continuous(breaks = seq(1970, 2020, 10)) +
  theme(plot.title = element_text(size =title.size, face = "bold", hjust=.5), 
        axis.text.x = element_text(size=x.axis.tick.size),  axis.title.x=element_blank(), 
        axis.text.y = element_text(size=y.axis.tick.size),  axis.title.y = element_text(size=y.axis.label.size),
        plot.margin=unit(c(0,3,0,3),"pt"))

#### Table 1 ####
cases_available <- unique(d.lit_attorney2$sample_id)

bal <- dockets_bj_summary1 %>%
  filter(year >= 2003) %>%
  mutate(in_sample = as.numeric(sample_id %in% cases_available))

attach(bal)
table(in_sample)
table(year, in_sample) %>%
  prop.table(margin = 1)
median(year[in_sample == 0]); median(year[in_sample == 1])
table(in_sample, cert_granted) %>%
  prop.table(margin = 1)
table(in_sample, ifp) %>%
  prop.table(margin = 1)
table(in_sample, doj_petitioners) %>%
  prop.table(margin = 1)
detach()

#### Table 2 ####
attach(d.lit_attorney2)
table(ifp)
table(ifp, cert_granted)
table(ifp, cert_granted) %>%
  prop.table(margin = 1)
table(ifp, pet_clerk)
table(pet_clerk, cert_granted, ifp)
table(pet_clerk, cert_granted, ifp) %>%
  prop.table(margin = c(1, 3))
detach()

#### Table 3 ####
#### Cert grant ~ former clerk on amicus briefs ####
# Crosstabs
attach(d.amc_docket1)
table(ifp)
table(ifp, cert_granted)
table(ifp, cert_granted) %>%
  prop.table(margin = 1)
table(ifp, brief_clerk_cat)
table(brief_clerk_cat, cert_granted, ifp)
table(brief_clerk_cat, cert_granted, ifp) %>%
  prop.table(margin = c(1, 3))
detach()

#### Table 4 ####
m.pet1 <- glm(cert_granted ~ pet_clerk,
              family = binomial(link = "logit"),
              d.lit_attorney2 %>%
                filter(ifp == 0))
rse1 <- coef_test(m.pet1, vcov = "CR1", 
                  cluster = d.lit_attorney2 %>%
                    filter(ifp == 0) %>%
                    pull(pet_attorney1),
                  test = "naive-t")$SE

m.pet2 <- glm(cert_granted ~ pet_clerk +
                doj_petitioners + opposition_respondents,
              family = binomial(link = "logit"),
              d.lit_attorney2 %>%
                filter(ifp == 0))
rse2 <- coef_test(m.pet2, vcov = "CR1", 
                  cluster = d.lit_attorney2 %>%
                    filter(ifp == 0) %>%
                    pull(pet_attorney1),
                  test = "naive-t")$SE

m.pet3 <- glm(cert_granted ~ pet_clerk +
                doj_petitioners + opposition_respondents +
                pet_lawschool_hy + pet_years_since_degree,
              family = binomial(link = "logit"),
              d.lit_attorney2 %>%
                filter(ifp == 0))
rse3 <- coef_test(m.pet3, vcov = "CR1", 
                  cluster = d.lit_attorney2 %>%
                    filter(ifp == 0) %>%
                    pull(pet_attorney1),
                  test = "naive-t")$SE

stargazer(m.pet1, m.pet2, m.pet3,
          se = list(rse1, rse2, rse3),
          dep.var.labels = "Cert granted",
          covariate.labels = c("Former clerk", "U.S. government", "Opposition brief",
                               "Harvard/Yale", "Career length"),
          star.cutoffs = c(0.05),
          star.char = c("*"),
          notes = "",
          notes.append = F,
          style = "ajps", type = "text")

#### Tables 5 ####
m.out <- matchit(pet_clerk ~ pet_lawschool + pet_years_since_degree_cat + doj_petitioners,
                 d.lit_attorney2 %>%
                   filter(ifp == 0),
                 method = "exact")
dm.lit_attorney <- match.data(m.out, group = "all")

# Pre-matching balance
bal_pet_attorney <- d.lit_attorney2 %>%
  filter(ifp == 0) %>%
  group_by(pet_clerk) %>%
  summarise(n_cases = n(),
            n_ifp = sum(ifp, na.rm = T),
            pct_ifp = mean(ifp, na.rm = T),
            n_doj = sum(doj_petitioners, na.rm = T),
            pct_doj = mean(doj_petitioners, na.rm = T),
            n_opp_filed = sum(opp_brief != "0_none", na.r = T),
            pct_opp_filed = mean(opp_brief != "0_none", na.rm = T),
            n_opp_clerk = sum(opp_brief == "2_clerk", na.rm = T),
            pct_opp_clerk = mean(opp_brief == "2_clerk", na.rm = T),
            n_lawschool_t14 = sum(pet_lawschool_t14, na.rm = T),
            pct_lawschool_t14 = mean(pet_lawschool_t14, na.rm = T),
            n_lawschool_hy = sum(pet_lawschool_hy, na.rm = T),
            pct_lawschool_hy = mean(pet_lawschool_hy, na.rm = T),
            q1_years_since_degree = quantile(pet_years_since_degree, .25),
            q2_years_since_degree = quantile(pet_years_since_degree, .5),
            q3_years_since_degree = quantile(pet_years_since_degree, .75)) %>%
  ungroup()

# Post-matching balance
bal_m.lit_attorney <- dm.lit_attorney %>%
  group_by(pet_clerk) %>%
  summarise(n_cases = n(),
            n_ifp = sum(ifp, na.rm = T),
            pct_ifp = mean(ifp, na.rm = T),
            n_doj = sum(doj_petitioners, na.rm = T),
            pct_doj = mean(doj_petitioners, na.rm = T),
            n_opp_filed = sum(opp_brief != "0_none", na.r = T),
            pct_opp_filed = mean(opp_brief != "0_none", na.rm = T),
            n_opp_clerk = sum(opp_brief == "2_clerk", na.rm = T),
            pct_opp_clerk = mean(opp_brief == "2_clerk", na.rm = T),
            n_lawschool_t14 = sum(pet_lawschool_t14, na.rm = T),
            pct_lawschool_t14 = mean(pet_lawschool_t14, na.rm = T),
            n_lawschool_hy = sum(pet_lawschool_hy, na.rm = T),
            pct_lawschool_hy = mean(pet_lawschool_hy, na.rm = T),
            q1_years_since_degree = quantile(pet_years_since_degree, .25),
            q2_years_since_degree = quantile(pet_years_since_degree, .5),
            q3_years_since_degree = quantile(pet_years_since_degree, .75)) %>%
  ungroup()

#### Table 6 ####
attach(dm.lit_attorney)
table(pet_clerk, cert_granted)
table(pet_clerk, cert_granted) %>%
  prop.table(margin = 1)
detach()

#### Table 7 ####
mm.pet1 <- glm(cert_granted ~ pet_clerk,
               family = binomial(link = "logit"),
               dm.lit_attorney)
rse1 <- coef_test(mm.pet1, vcov = "CR1", 
                  cluster = dm.lit_attorney$pet_attorney1, test = "naive-t")$SE

mm.pet2 <- glm(cert_granted ~ pet_clerk + pet_lawschool + pet_years_since_degree,
               family = binomial(link = "logit"),
               dm.lit_attorney)
rse2 <- coef_test(mm.pet2, vcov = "CR1", 
                  cluster = dm.lit_attorney$pet_attorney1, test = "naive-t")$SE

stargazer(mm.pet1, mm.pet2,
          se = list(rse1, rse2),
          dep.var.labels = "Cert granted",
          omit = "pet_lawschool.*",
          covariate.labels = c("Former clerk", "Career length"),
          star.cutoffs = c(0.05),
          star.char = c("*"),
          notes = "",
          notes.append = F,
          style = "ajps", type = "text")

#### Figure 3 ####
attach(dm.lit_attorney)
newdata <- tibble(pet_lawschool_hy = 1,
                  pet_lawschool = "Harvard",
                  pet_years_since_degree = mean(pet_years_since_degree))
detach()

pred <- bind_rows(
  ggpredict(mm.pet1, "pet_clerk",
            vcov.fun = "vcovCR", vcov.type = "CR1", 
            vcov.args = list(cluster = dm.lit_attorney$pet_attorney1)) %>%
    tibble() %>%
    mutate(model = "Model 1"),
  ggpredict(mm.pet2, "pet_clerk",
            condition = c(pet_lawschool = "Harvard",
                          pet_years_since_degree = newdata$pet_years_since_degree)) %>%
    tibble() %>%
    mutate(model = "Model 2")
) %>%
  rename(attorney_type = x)

ggplot(pred, aes(x = model, y = predicted, color = factor(attorney_type))) + theme_bw() +
  geom_point(size=5) +   scale_color_manual(values = c("gray", "black")) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, size = 1) +
  annotate("text", x = 1.5, y = .16, label = "Former clerks", size=5) +
  annotate("text", x = 1.5, y = .06, label = "Non-former clerks", size=5) +
  labs(x = "", y = "Predicted Probability of Cert Grant") +
  coord_cartesian(ylim=c(0,.25)) +  scale_y_continuous(expand = c(0, 0)) +
  theme(plot.title = element_text(size =title.size, face = "bold", hjust=.5), 
        axis.text.x = element_text(size=x.axis.tick.size),  axis.title.x=element_text(size=y.axis.label.size), 
        axis.text.y = element_text(size=y.axis.tick.size),  axis.title.y = element_text(size=y.axis.label.size), legend.position = "none")

#### Table 8 ####
# Simple logit
m.amc1 <- glm(cert_granted ~ brief_clerk_cat,
              family = binomial(link = "logit"),
              d.amc_docket1 %>%
                filter(ifp == 0))
linearHypothesis(m.amc1, c("brief_clerk_cat1_brief_with_clerk=brief_clerk_cat2_brief_without_clerk"), test="F")

m.amc2 <- glm(cert_granted ~ brief_clerk_cat +
                doj_petitioners + pet_clerk,
              family = binomial(link = "logit"),
              d.amc_docket1 %>%
                filter(ifp == 0))
linearHypothesis(m.amc2, c("brief_clerk_cat1_brief_with_clerk=brief_clerk_cat2_brief_without_clerk"), test="F")

m.amc3 <- glm(cert_granted ~ n_amc_clerks + n_amc_nonclerks +
                n_briefs_for + doj_petitioners + pet_clerk,
              family = binomial(link = "logit"),
              d.amc_docket1 %>%
                filter(ifp == 0))
linearHypothesis(m.amc3, c("n_amc_clerks = n_amc_nonclerks"), test="F")

m.amc4 <- glm(cert_granted ~ n_clerks_per_brief + n_nonclerks_per_brief +
                n_briefs_for + doj_petitioners + pet_clerk,
              family = binomial(link = "logit"),
              d.amc_docket1 %>%
                filter(ifp == 0))
linearHypothesis(m.amc4, c("n_clerks_per_brief = n_nonclerks_per_brief"), test="F")

stargazer(m.amc1, m.amc2, m.amc3, m.amc4,
          dep.var.labels = "Cert granted",
          covariate.labels = c("Amicus brief with former clerks", "Amicus brief without former clerk",
                               "Number of former clerks on brief", "Number of non-former clerks on brief",
                               "Number of former clerks per amicus brief", "Number of non-former clerks per amicus brief",
                               "Number of amicus briefs",
                               "U.S. government", "Petition by former clerk"),
          star.cutoffs = c(0.05),
          star.char = c("*"),
          notes = "",
          notes.append = F,
          style = "ajps", type = "text")

#### Figure 4 ####
attach(d.amc_docket1)
newdata <- tibble(ifp = 0,
                  doj_petitioners = table(doj_petitioners) %>%
                    sort(decreasing = T) %>% names() %>% .[1] %>% as.numeric(),
                  pet_clerk = table(pet_clerk) %>%
                    sort(decreasing = T) %>% names() %>% .[1] %>% as.numeric())
detach()

pred <- bind_rows(
  ggpredict(m.amc1, "brief_clerk_cat") %>%
    tibble() %>%
    mutate(model = "Model 1"),
  ggpredict(m.amc2, "brief_clerk_cat",
            condition = c(ifp = newdata$ifp,
                          doj_petitioners = newdata$doj_petitioners,
                          pet_clerk = newdata$pet_clerk)) %>%
    tibble() %>%
    mutate(model = "Model 2")
) %>%
  rename(brief_type = x)

ggplot(pred,
       aes(x = model, y = predicted, color = factor(brief_type))) +
  geom_point(size = 3) +
  scale_color_manual(values = c("gray80", "gray50", "black")) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
                width = 0, size = 1) +
  annotate("text", x = 1.5, y = .03, label = "No brief") +
  annotate("text", x = 1.5, y = .25, label = "Without former clerks", angle=-10) +
  annotate("text", x = 1.5, y = .45, label = "With former clerks", angle=-10) +
  labs(x = "", y = "Predicted Probability of Cert Grant") + 
  coord_cartesian(ylim=c(0,.65)) +  scale_y_continuous(breaks=seq(0,.6,.1), expand = c(0, 0)) +
  theme_bw() + 
  theme(plot.title = element_text(size =title.size, face = "bold", hjust=.5), 
        axis.text.x = element_text(size=x.axis.tick.size),  axis.title.x=element_text(size=y.axis.label.size), 
        axis.text.y = element_text(size=y.axis.tick.size),  axis.title.y = element_text(size=15), legend.position = "none")

#### Table 9 ####
attach(d.lit_attorney2 %>%
         filter(ifp == 0,
                doj_petitioners == 0,
                doj_respondents == 0))
table(pet_clerk, opp_brief)
table(pet_clerk, cert_granted, opp_brief)
table(pet_clerk, cert_granted, opp_brief) %>%
  prop.table(margin = c(1, 3))
detach()

#### Table A-1 ####
attach(d.lit_attorney2 %>%
         filter(pet_clerk == 1, ifp == 0))
table(pet_sitting_justice, cert_granted)
table(pet_sitting_justice, cert_granted) %>%
  prop.table(margin = 1)
detach()

#### Table A-2 ####
m.pet4 <- glm(cert_granted ~ pet_sitting_justice,
              family = binomial(link = "logit"),
              d.lit_attorney2 %>%
                filter(pet_clerk == 1,
                       ifp == 0))
rse4 <- coef_test(m.pet4, vcov = "CR1", 
                  cluster = d.lit_attorney2 %>%
                    filter(pet_clerk == 1,
                           ifp == 0) %>%
                    pull(pet_attorney1),
                  test = "naive-t")$SE

m.pet5 <- glm(cert_granted ~ pet_sitting_justice + doj_petitioners +
                pet_lawschool_hy + pet_years_since_degree_cat,
              family = binomial(link = "logit"),
              d.lit_attorney2 %>%
                filter(pet_clerk == 1,
                       ifp == 0))
rse5 <- coef_test(m.pet5, vcov = "CR1", 
                  cluster = d.lit_attorney2 %>%
                    filter(pet_clerk == 1,
                           ifp == 0) %>%
                    pull(pet_attorney1),
                  test = "naive-t")$SE

stargazer(m.pet4, m.pet5,
          se = list(rse4, rse5),
          dep.var.labels = "Cert granted",
          covariate.labels = c("Sitting Justice", "U.S. government", "Harvard/Yale", "Career length"),
          star.cutoffs = c(0.05),
          star.char = c("*"),
          notes = "",
          notes.append = F,
          style = "ajps", type = "text")

#### Table A-3 ####
# Cert grant rates within non-missing and missing data
attach(d.lit_attorney_add)
table(ifp, cert_granted, case_available)
table(ifp, cert_granted, case_available) %>%
  prop.table(margin = c(1, 3))
detach()

# Simple logit including sample of missing data (attorney background)
m.pet1_add <- glm(cert_granted ~ pet_clerk,
                  family = binomial(link = "logit"),
                  d.lit_attorney_add %>%
                    filter(ifp == 0))
rse1_add <- coef_test(m.pet1_add, vcov = "CR1", 
                      cluster = d.lit_attorney_add %>%
                        filter(ifp == 0) %>%
                        pull(pet_attorney1),
                      test = "naive-t")$SE

m.pet2_add <- glm(cert_granted ~ pet_clerk +
                    doj_petitioners + opposition_respondents,
                  family = binomial(link = "logit"),
                  d.lit_attorney_add %>%
                    filter(ifp == 0))
rse2_add <- coef_test(m.pet2_add, vcov = "CR1", 
                      cluster = d.lit_attorney_add %>%
                        filter(ifp == 0) %>%
                        pull(pet_attorney1),
                      test = "naive-t")$SE

m.pet3_add <- glm(cert_granted ~ pet_clerk +
                    doj_petitioners + opposition_respondents +
                    pet_lawschool_hy + pet_years_since_degree,
                  family = binomial(link = "logit"),
                  d.lit_attorney_add %>%
                    filter(ifp == 0))
rse3_add <- coef_test(m.pet3_add, vcov = "CR1", 
                      cluster = d.lit_attorney_add %>%
                        filter(ifp == 0) %>%
                        pull(pet_attorney1),
                      test = "naive-t")$SE

stargazer(m.pet1_add, m.pet2_add, m.pet3_add,
          se = list(rse1_add, rse2_add, rse3_add),
          dep.var.labels = "Cert granted",
          covariate.labels = c("Former clerk", "U.S. government", "Opposition brief",
                               "Harvard/Yale", "Career length"),
          star.cutoffs = c(0.05),
          star.char = c("*"),
          notes = "",
          notes.append = F,
          style = "ajps", type = "text")

#### Table A-4 ####
m.out <- matchit(pet_clerk ~ pet_lawschool + pet_years_since_degree_cat + doj_petitioners,
                 d.lit_attorney_add %>%
                   filter(!is.na(pet_lawschool),
                          !is.na(pet_years_since_degree_cat),
                          ifp == 0),
                 method = "exact")
dm.lit_attorney_add <- match.data(m.out, group = "all")

mm.pet_add1 <- glm(cert_granted ~ pet_clerk,
                   family = binomial(link = "logit"),
                   dm.lit_attorney_add)
rse1 <- coef_test(mm.pet_add1, vcov = "CR1", 
                  cluster = dm.lit_attorney_add$pet_attorney1, test = "naive-t")$SE

mm.pet_add2 <- glm(cert_granted ~ pet_clerk + pet_lawschool_hy + pet_years_since_degree,
                   family = binomial(link = "logit"),
                   dm.lit_attorney_add)
rse2 <- coef_test(mm.pet_add2, vcov = "CR1", 
                  cluster = dm.lit_attorney_add$pet_attorney1, test = "naive-t")$SE

stargazer(mm.pet_add1, mm.pet_add2,
          se = list(rse1, rse2),
          dep.var.labels = "Cert granted",
          covariate.labels = c("Former clerk", "Harvard/Yale", "Career length"),
          star.cutoffs = c(0.05),
          star.char = c("*"),
          notes = "",
          notes.append = F,
          style = "ajps", type = "text")